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We consider discs that orbit a central object and are tidally perturbed by a circular 
orbit companion. Such discs are sometimes subject to an eccentric instability due to the 
effects of certain resonances. Eccentric instabilities may be present in planetary rings 
perturbed by satellites, protostellar discs perturbed by planets, and discs in binary star 
systems. Although the basic mechanism for eccentric instability is well understood, 
the detailed response of a gaseous disc to such an instability is not understood. We 
apply a linear eccentricity evolution equation developed by Goodchild and Ogilvie. We 
explore how the eccentricity is distributed in such a disc and how the distribution in 
turn affects the instability growth rate for a range of disc properties. We identify a disc 
mode, termed the superhump mode, that is likely at work in the superhump binary 
star case. The mode results from the excitation of the fundamental free precession 
mode. We determine an analytic expression for the fundamental free mode precession 
rate that is applicable to a sufficiently cool disc. Depending on the disc sound speed 
and disc edge location, other eccentric modes can grow faster than the superhump 
mode and dominate. 

Key words: accretion, accretion discs — instabilities — (stars:) binaries (including 
multiple): close — planet disc interactions — planets: rings — methods: numerical 



Consider a gaseous disc that orbits about a central object. If 
the disc is tidally perturbed by a circular orbit companion, 
but the disc is otherwise circular, then its response at certain 
(Lindblad) resonance locations can be quite strong and re- 
sult in the launching of waves that extract energy and angu- 
lar momentum from the companion (Goldreich & Tremaine 
1979). If the perturber were to be removed, the disc would 
settle back to a circular state, with its mass having been 
rearranged from its initial state by the resonant torques. 
For discs that are slightly eccentric, additional resonances 
arise. Some of these resonances have the property that they 
can cause the disc eccentricity to grow exponentially in time, 
even though the perturber is on a circular orbit. In this case, 
if the perturber were to be removed, the disc would be in 
an eccentric state. The process by which eccentricity growth 
occurs in a fluid disc can be understood through a mode- 
coupling mechanism (Lubow 1991a; hereafter L91). The in- 
stability can also be understood in terms of particle dy- 
namics with dissipation (Borderies, Goldreich, & Tremaine 
1983, Ogilvie 2007). Such instabilities may explain the ec- 
centricities of planetary rings perturbed by satellites (Bor- 
deries et al 1983), the eccentricities of circumstellar discs 



perturbed by planets or brown dwarfs (Papaloizou, Nelson, 
& Masset 2001, Kley & Dirksen 2005, D'Angelo, Lubow, & 
Bate 2006), and eccentricities in discs of catalysmic binaries 
(Whitehurst 1988, Lubow 1991b, Osaki 1996), and X-ray 
binaries (Haswell et al. 2001, Neil, Balyn, & Cobb 2007). 

In this paper, we will mainly consider the superhump 
binary case which has been studied, both observationally 
and theoretically. We will generalise the results somewhat 
to cover cases involving warmer (larger dimensionless disc 
thickness) protostellar discs. The initial observational evi- 
dence for eccentric discs came from an analysis of a drifting 
feature, the superhump, in the hght curves of extreme mass 
ratio cataclysmic binaries that undergo unusually strong 
outbursts called superoutbursts (Warner 1975). The period 
of the superhump feature is slightly longer than the binary 
orbit period. The drifting feature in the light curve was at- 
tributed to an eccentric prograde precessing disc (Vogt 1982, 
Osaki 1985). The prograde precession is due to the grav- 
itational effects of the companion that causes departures 
of disc streamlines from closed Keplerian ellipses. Weaker 
retrograde effects occur due to gas pressure (Lubow 1992, 
Murray 2000). Simulations by Whitehurst (1988) revealed 
an eccentric instability that can be understood through the 
effects of the 3:1 eccentric Lindblad resonance (L91). Many 
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observational studies liave supported the idea that the disc 
is eccentric (e.g., Osaki 2003, Patterson 2005). 

Several particle based simulations (SPH or something 
similar) of discs in binaries revealed properties of the insta- 
bility and the disc precession (see, e.g.. Smith et al 2007). 
For many years, simulations carried out by grid-based codes 
were unable to find this instability and cast some doubt 
on the eccentric disc model (e.g., Heemskerk 1994, Stehle 
1999). However, Kley, Papaloizou, & Ogilvie (2008) recently 
reported finding this instability with a grid-based code and 
showed that its properties largely agree with theoretical 
expectations of the mode-coupling model. However, some 
properties of the instability revealed by these simulations 
were not predicted by the L91 model, such as the depen- 
dence of the growth rate on the disc sound speed. 

The L91 model determined the eccentricity growth rate 
for a ring of gas having uniform eccentricity. It was gener- 
alised to a full disc, again assuming the eccentricity is uni- 
formly distributed. However, it is by no means obvious that 
the disc eccentricity would always be uniform. The distribu- 
tion of eccentricity has an important influence on the eccen- 
tricity growth rate. The resonance provides a certain rate of 
angular momentum and energy input, resulting in an eccen- 
tricity growth. But if this angular momentum and energy 
input is shared over a broad eccentric region of the disc, 
then we expect the eccentricity growth rate to be smaller 
than if the region is narrow. In fact, the L91 model predicts 
that the growth rate to vary as roughly the inverse of the 
radial width of the eccentric region. 

Important progress in understanding the disc eccentric- 
ity distribution and its effects on the growth rate was made 
by Goodchild & Ogilvie (2006; hereafter GO06). They de- 
rived differential equations for the radial distribution and 
growth rate of eccentricity, subject to the eccentric insta- 
bility caused by the resonance. Their results for superhump 
binaries, involving the 3:1 resonance, indicated that the disc 
eccentricity growth rate was lower than expected from sim- 
ulations and observations. The reason for the small growth 
rate is that the mode they analysed has a low amplitude near 
the resonance. Consequently, the resonance was ineffective 
in providing eccentricity growth. We therefore have a puzzle 
that a more complete modeling of the disc eccentricity evo- 
lution resulted in growth rates that are lower than expected. 
In this paper we re-examine the nature of the eccentricity 
evolution of a disc, using the GO06 equations, in an attempt 
to resolve this issue. In the process, we determine a more 
strongly excited mode that we call the superhump mode. 
We analyse its properties under various conditions. We also 
determine other modes than can dominate, under somewhat 
different disc conditions than are typically expected in su- 
perhump binaries. 

The outline of the paper is as follows. In Section 2, we 
describe the eccentricity evolution equations. Section 3 spec- 
ifies the standard model parameters used in the calculations. 
Section 4 discusses the methods to solve the eccentricity evo- 
lution equations. The results of the calculations for a variety 
of system parameters are described in Section 5. Section 6 
contains the summary. 



2 ECCENTRICITY EVOLUTION EQUATIONS 

Consider a two-dimensional disc that orbits about a primary 
object which is located at the origin of a cylindrical coordi- 
nate system (r, 4>) . The low-mass secondary is in a circular 
orbit of semi-major axis d about the central object. We ap- 
ply the linearised eccentricity equations of GO06, expressed 
in the following form 

idr{air)drE{r,t)) + ib{r)E{r,t) + J{r)s{r)E{r,t) 

= J{r)dtEir,t), (1) 

where E{r,t) — e{r,t) exp {izAj{r,t)) is the complex eccen- 
tricity, for real eccentricity e and periapse angle zu. This two 
dimensional approximation is valid under the Eissumption of 
vertical (perpendicular to the disc plane) hydrostatic equi- 
librium for a thin disc. The reduction to two dimension re- 
sults from a vertically integration of three dimensional equa- 
tions. We ignore the effects of viscosity on the eccentricity. 
For the case of no perturbing body, the nonlinear evolution 
of an eccentric, three dimensional, viscous disc was analysed 
by Ogilvie (2001). 

E is related to the linear perturbations from the axisym- 
metric circular velocity. For the velocity expressed cylindri- 
cal coordinates as {u {r,t) exp {—i(f)),v' {r,t) exp {—ifj))) we 
have that 

u'(r,t) = irn{r)E{r,t) (2) 

and 

v'{r,t) = ^rQir)E{r,t), (3) 

where 0(r) is the Keplerian orbital frequency about the pri- 
mary of mass Mp given by 

Quantity J(r) is the disc angular momentum per unit radius 
divided by tt and is given by 

J(r) = 2-r''0('r)S(r). (5) 

Functions a{r) and 6(r) are given by 

a{r) = 7P(r)^^ (6) 

b{r) = ^r' + J{r)^„ (7) 
ar 

where P{r) is the two-dimensional (vertically integrated) 
disc pressure, S(r) is the disc surface density, and -djg is the 
gravitational precession rate of a free particle on an eccentric 
orbit which is given by 

^^^'^^ = 1^^ ^3)2 O- (8) 

The terms involving quantities a and b describe the eccen- 
tricity propagation and precession, respectively. Quantity 7 
is the gas adiabatic index. Mass ratio q is the mass of the 
secondary (perturbing) object divided by the mass of the 
primary (central) object, and feg^j is the Laplace coefficient 
for the m — 1 tidal potential component associated with 
the perturbing object. Real function s(r) is the eccentricity 
growth rate contribution from the eccentric instability. It is 
largest in a region containing the resonance that drives the 
instability. In the case of superhump binaries, this is the 3:1 
resonance. 
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The surface density E(r-) is taken to vary smoothly from 
inner radius ri to outer radius ro and to be zero outside 
this range. Functions a{r) and b{r) vary smoothly in radius, 
while function s(r) is localised near the resonance. For the 
eccentricity injection rate caused by a resonance, we gener- 
ally adopt a gaussian form 

X exp [-(r - r^osf/wl 



s(r) = 



2 1 

'^rcsj 



(9) 



where rrcs is the radius of the resonance that excites eccen- 
tricity, Wres is the resonance width, and x is a- measure of 
the resonance strength. 

We consider the boundary conditions at the disc inner 
and outer edges located at radii ri and ro, respectively. Fol- 
lowing GO06, we adopt the outer boundary condition 



drE{ro,t) = 0. 



(10) 



For a Keplerian disc, the divergence of the velocity is pro- 
portional to drE. Consequently, this boundary condition is 
equivalent to requiring that the Lagrangian density pertur- 
bation near the disc outer edge vanishes. 

The inner boundary condition is taken to be 



S(n,t) = 0. 



(11) 



This condition could be imposed by a hard wall inner circu- 
lar boundary due to a central object. We discuss the inner 
boundary condition further in Section [5.71 



3 STANDARD MODEL 

We start with a simple disc model that we define here. We 
later explore certain variations from this model. We adopt 
disc models in which 7 = 5/3, T oc r~^^*, and E oc r~^^^. 
The density form corresponds to a steady state alpha disc 
with a constant value of alpha. In this model, the disc 
thickness-to-radius ratio is given by 

H/r = h{r/r,,,y^\ (12) 

where h — H/r evaluated at r = rrcs. We do not consider 
any tapering of the disc density near the the outer edge. 
Some tapering is needed over a distance ~ H to prevent 
pressure forces from causing the disc to become Rayleigh 
unstable, i.e., in order that the epicyclic frequency be real. 
We do not include tapering, in order to keep the standard 
model as simple as possible. In this paper, we consider discs 
with various h values, but take the r— dependence of T and 
E to be fixed. 

We take the unit of length to be the binary separation 
d. We consider the effects the eccentricity instability at the 
3:1 eccentric Lindblad resonance. For q small, the resonant 
radius is approximately given by 

r.es = 3-'/'(l + g)-^/^ (13) 

The 3:1 resonance strength x in equation ((9]) is given by 



X = 2.08(7 ^^brros, 



(14) 



where ilb is the binary orbital frequency (L91). This expres- 
sion was independently confirmed by the particle model of 
Ogilvie (2007). The resonance width Wros depends on the 
disc sound speed and we take 

Wrcs — tl'^^ r-rea (15) 



(see Meyer- Vernet & Sicardy 1987). We consider variations 
from this value of the width in Section 15.81 

We adopt a mass ratio g = 0.1 as a standard value. With 
this mass ratio, the disc outer radius determined by Paczyn- 
ski's (1977) particle orbit intersection method is equal to 
0.46. The disc radius given by the prescription suggested 
in Whitehurst & King (1991) evaluates to 0.53, which is 
0.9 times the Roche lobe radius. We adopt a fiducial value 
of 0.5. The disc inner radius is chosen as 0.01 for numeri- 
cal convenience to avoid the singularity at the disc center. 
Smaller inner disc radius values have little impact on the 
results presented here, as is discussed in Section [5.71 

In summary, we adopt a set of standard parameter val- 
ues to model superhump binaries and later consider varia- 
tions from some of these values. The standard model has 
q = 0.1, h = 0.02, n = 0.01, ro = 0.5, r^os = 0.466, 7 = 5/3, 
T oc r-3''^ and E oc r~^/*. 



4 METHODS OF SOLUTION 

Equation ([T]) and boundary conditions (|10[) and (|11[) are 
homogeneous and therefore permit the specification of the 
value of E{r, 0) at some selected radius r. For the standard 
model we adopt a normalization 



S(ro,0) 



(16) 



Of course, we are not suggesting that the actual eccentricity 
is unity. But, since the equations are linear, the solutions 
can be scaled linearly to any desired normalization. 

Equation ([T]) can be solved by searching for eigenmodes 
in which 



E{r,t) = E{r)exp (iujt) 



(17) 



and complex frequency cj is the eigenvalue. The function 
E{r) and eigenfrequency uj can be determined by shooting 
methods that satisfy the boundary and normalization con- 
ditions. Equation ^ is second order in space and can be 
integrated inward in radius starting at r = ro by using the 
two conditions, equations (jlOp and H16[) . and an assumed 
value for lj. Complex frequency uj is adjusted until the in- 
ner boundary condition, equation (|11|) . is satisfied by using 
Newton's method. Another approach would be to discretise 
equation ^ in radius and solve the equations represented 
by a sparse matrix. 

Such approaches have the drawback that there are many 
modes in the system and one must search for the physi- 
cally relevant one that is the fastest growing. The shooting 
method is particularly prone to missing the fastest growing 
mode. The reason is that the guessed value for the eigen- 
frequency must be close to the value for the fastest growing 
mode. Otherwise, the method will typically converge on an- 
other eigenfrequency of the system. 

We have instead opted to find the fastest growing mode 
by integrating equation ([TJ in space and time as a PDE. We 
apply the following initial condition 



E{r,0) = ( 



(18) 



for ri ^ r ^ To. This initial eccentricity satisfies boundary 
and normalization conditions (|10p . (Illf) . and (|16l) . Over time 
E{r,t)/ E(ro,t) settles to an eigenfunction in which dtE/E 
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is constant in space and is the complex eigenfrequency. This 
approach is computationally much slower than solving the 
eigenvalue problem with the shooting method, if one has 
a good initial guess for the eigenfrequency. We apply the 
eigenfrequency obtained from the time integration to the 
shooting method described above to verify its existence. We 
then smoothly vary other parameters such as q and h to 
rapidly obtain many other solutions by iteratively apply- 
ing the eigenfrequency from the last obtained solution as 
the initial guess of the eigenfrequency for the next nearby 
set of parameters. In this way, we obtain a set of solutions, 
along a continuous branch that began with a fastest growing 
mode. There is no guarantee, however, that all solutions in 
this branch are the fastest growing for their set of parame- 
ters. Sufficiently nearby solutions are typically found to also 
be the fastest growing. However, as we will see, multiple 
branches occur that contain the fastest growing modes in 
different parameter regimes. 

To integrate equation ([ij, we initially wrote a C pro- 
gram to implement the Crank-Nicholson method (e.g., 
Watanabe & Tsukada 2000). This method is unconditionally 
stable and fast, but produced some unwanted oscillations 
at high resolution, possibly as a consequence of Gudunov's 
Theorem (Wesseling 2001). Instead, we opted for the method 
of lines in Mathematica. With this method, the eccentric- 
ity is discretised in radius, producing a set of eccentricities 
Ei{t), one at each radial grid point i. Equation ((T)) becomes 
to a set of coupled ODEs in time that are integrated by 
high order schemes. This method produced high resolution 
results. 

Based on G096, we express the growth rate for the 
standard model as 

J r^/'^e'^{r)dr max{we,Wrcs) \ e^ax J 

(19) 

where the integrals are taken over the disc from r — ri to 
r = To, emax is the maximum eccentricity in the disc, and 
We is the radial width of the eccentricity distribution. This 
equation demonstrates the dependence of the growth rate 
on the properties of the eccentricity distribution. 



5 RESULTS 

5.1 Superhump Mode 

The method of lines solution to equation lU for the stan- 
dard model, subject to initial condition (|18p and boundary 
conditions (|10|l and (|lip . is plotted in Fig. [T] After a time 
of about 16 binary orbital periods, the eccentricity distri- 
bution settles into a mode. The values of the normalised 
eccentricity e{r,t)/e{ro,t) and the phase angle ■ou{r,t) be- 
come nearly time-independent. The eccentricity growth rate 
dte{r,t) /e{r,t) and precession rate dtvj{r,t) become inde- 
pendent of radius. 

The eccentricity growth rate for the standard model is 
determined to be 

dtey,t) ^ ^j^(^) ^ o^oi^ (20) 
e{r,t) ^ ' ^ ' 

where cj = dtE{r,t)/ E{r,t). 

This rapid growth rate is a consequence of the strong 



overlap between the eccentricity distribution in the mode 
with the eccentricity driving by the resonance and the nar- 
rowness of the eccentricity distribution. From the figure, we 
estimate We — 0.1 and e(rres) — Smax- Applying the param- 
eters for the standard model {q — 0.1 and rros — 0.466), 
we obtain an estimated growth rate from equation (|19p that 
is —Im{u)) ~ 0.1f2b, in agreement with the rate obtained 
from the time-evolution solution, equation (|20|l . Therefore, 
we find that the eccentricity is robustly excited in the outer 
parts of the disc. This calculation is telling us something 
that L91 did not, the effects of the eccentricity distribution 
on the growth rate. As we will see, the width of the eccen- 
tricity distribution is a function of the system parameters 
and the growth rate itself. 

The mode shown in the figure is termed the super- 
hump mode. To understand its properties, we associate a 
timescale with each coefficient on the left-hand side of equa- 
tion IT}. The pressure, precession, and resonance instability 
timescales are estimated as vj^J/a, J/b and fo/x, respec- 
tively. For the standard model, they crudely evaluate near 
the resonance to IG^ vol / {r^Q.i,) , '20/O,b, and 50/Q,b, respec- 
tively. The effects of pressure act to spread the eccentricity 
distribution and balance against the effects of precession and 
resonance instability that will be shown to confine the ec- 
centricity. The eccentricity distribution width We must be 
substantially smaller than ro, so that the pressure timescale 
is comparable to the precession and instability timescales. 
Therefore, the eccentricity should be radially confined to a 
region smaller than the disc radius. 

Notice that the phase of the mode decreases with radius, 
indicating that the mode is a trailing waveform. The radial 
wavelength is long, about equal to 0.3. Recall that the disc 
scale height is only about 2% of the disc radius or 0.01. 

As outlined in Section|4l we apply the complex eigenfre- 
quency obtained for the standard model to solve the mode 
equation by the shooting method. We then smoothly vary 
parameters to determine how the growth rate and modal 
structure varies as a function of various parameters, as will 
be described in the subsequent subsections. 

5.2 Variations in q 

We consider how the growth rate varies with binary mass 
ratio, q. As q varies, the radius of the 3:1 resonance varies 
slightly, according to equation (|15p . But since the resonance 
lies close to the disc outer edge, varying q alone would intro- 
duce an additional effect in changing the clearance between 
resonance and disc edge that we denote as 

Arc = ro — rros. (21) 

We analyse the effects of changing the clearance Arc sepa- 
rately. Consequently, when varying q we also vary ro so as to 
maintain a constant value Arc = 0.5 — 0.466 = 0.034. This 
means that ro varies slightly with q. 

We carried out a sequence of solutions to the modal 
equations using the shooting method for different binary 
mass ratios q. The results plotted in Fig. [2] show that the 
growth rate varies somewhat faster than q^. A quadratic de- 
pendence of the growth rate on q follows from the x term 
(see equation (|14|l ) in equation (|19p . The somewhat more 
rapid dependence is a consequence of additional mode con- 
finement with increasing q (see Fig. |3]). The confinement is 
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Figure 1. Evolution of eccentricity for tlie standard model (q = 0.1,h = 0.02, To = 0.5). The initial eccentricity distribution is of a cosine 
form, equation ||18|I . Radius r is in units of the binary separation. Left panel: e(r) at various times, normalised by its value at the disc 
outer edge. The curves montonically drop in time at r = 0.2, corresponding to times of 0, 1, 2, 8 and 16 (thick line) binary orbits. The 
dotted curve is the 3:1 resonance instability factor s(r) of equation l(lj in units of Ob- Right panel: Phase angle (periapse angle) zuir) of 
the eccentricity in radians divided by tt as a function of radius at various times. The curves monotonically increase in time at r = 0.4, 
corresponding to times of (fiat line of zero phase), 2, 8, and 16 (thick line) binary orbits. 




Figure 3. Eccentricity distributions for cases with standard pa- 
rameters, but having different values of binary mass ratio q. The 
q values from the highest to lowest solid curves are 0.05, 0.1, and 
0.2, respectively. The disc outer edge varies slightly with q from 
the standard value of 0.5, so that the distance from the resonance 
to the disc outer edge is maintained at a constant value. The dot- 
ted curves are the 3:1 resonance instability factors s(r) of equation 
JTJ in units of O]-, for the same set of q values, where the peak 
values increase with increasing q. The eccentricity distributions 
narrow with increasing q, since the increasing growth rates cause 
the eccentricity to be more concentrated near the resonance. 



a result of a competition between the pressure that spreads 
the distribution with the eccentricity growth that narrows 
the distribution to a region near the resonance. 



5.3 Free eccentric mode 

To determine the free eccentric mode corresponding to the 
superhump mode, we regarded q as fixed and constructed a 
sequence of models starting with the superhump mode and 
smoothly reduced the resonance strength x (or s(r)) to zero. 
In that limit, we obtain a mode structure that is plotted 
in Fig. U along with the mode structure in the standard 



model. Both modes are computed with the same disc outer 
radius equal to 0.5. The x ~ ^ mode corresponds to the 
fundamental (lowest order) free eccentric mode. Therefore, 
we see that the superhump mode results from the excitation 
of the free eccentric mode. 

The additional mode confinement in the case of the 
growing mode can again be understood in terms of a com- 
petition between the radial pressure-induced propagation of 
eccentricity that acts to spread the mode and its growth 
that acts to localise the mode close to the 3:1 resonance 
that resides near the disc outer edge (the dotted curve in 
the figure). The confinement in the free mode case is shown 
below to be due to the effects of disc precession. 

A previous estimate of eccentric disc free precession 
rate, based on the WKB theory of density waves, suggested 

a; = t^g(ro)- ^r!(ro), (22) 

where k is the radial wavenumber of the eccentric mode 
(Lubow 1992). This result demonstrates that pressure causes 
retrograde precession (by the second term on the right-hand 
side), opposite to the prograde precession due to the com- 
panion (the first term on the right-hand side). However, the 
magnitude of k was not determined. 

We analyse the fundamental free precession mode for a 
cool disc by means of equation together with equation 
(|17|l and the condition s(r) = 0. For a cool disc, we expect 
that the mode is confined to a small radial region near the 
disc outer edge and so \dE/dr\ » \E\/r. In addition, we 
expect that uj ~ ri7g(r-o). We apply these approximations to 
obtain 

(P E 

a(r-o)-^ + ■^(^o)(^g('") - '-')E{r) = 0. (23) 

We expand ri7g(r) in a Taylor series to linear order about To 
and then have 

(24) 



where Ar — r - 



Co 



(25) 
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Figure 2. Left panel: Eccentricity growth rate in units of Q^,. plotted against binary mass ratio q for the standard model. The disc outer 
edge varies slightly with q from the standard value of 0.5, so that the distance from the resonance to the disc outer edge is maintained 
at a constant value. Right panel: Plot of 13, defined by din (e)/dt oc as a function of q based on the results of the left panel. 
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Figure 4. Eccentricity distribution for a system with standard 
parameters (solid line) and a system with zero resonance strength 
(dashed line) , obtained through a sequence of models in which the 
resonance strength of the standard parameter case is smoothly 
lowered to zero. The dotted line is the 3:1 resonance instability 
factor sir) of equation JlJ in units of Ob for the standard case 
(solid line). The 3:1 resonance in the outer parts of the disc acts 
to further confine the mode near the disc outer edge. The dashed 
curve is the fundamental free eccentric mode of the disc. It is that 
mode which is excited by the eccentricity instability to become 
the superhump mode. 



in order that the outer boundary condition (|10p is satisfied. 
There are many (formally infinitely many) values of C2 that 
satisfy equation (|29[1. since Ai is oscillatory. Each such so- 
lution corresponds to a free eccentric mode. The number 
of nodes increases with increasing values of C2. The funda- 
mental mode has the smallest value of C2 and for this mode 
C2 « 1.019. 

The precession rate is given by 

a; = ri7g(ro) — C2 lUe (30) 
ar 

or 

a; ~ •ri7g(ro — C2Uie). (31) 

Therefore, the precession rate for a cool disc with pressure 
is equal to the gravitational precession rate evaluated at a 
distance C2We — We inside the disc edge. 

Equation (|28p for width We demonstrates the competi- 
tion between the sound speed or H that broadens the eccen- 
tricity distribution, and differential gravitational precession, 
dtjo^jdr, that acts to confine it. The analytic eccentricity 
distribution given by equation (|27|) roughly agrees with the 
numerically determined distribution for Hjr = 0.02, but be- 
comes fairly accurate for Hjr — 0.005 (see Fig. [5]). The pre- 
cession rates given by equation (|30} obtain similar levels of 
agreement with numerically determined values (see Fig. 



and 



2 / r \ 2 1 d-ri7„ 



Cl = 



-y \ HJ n dr 



(26) 



where Ci and C2 are evaluated at r = Tq. Applying equations 
(|11|) and (|16|) . we obtain the solution to equation (|24|) as 



Ai{-C2 - (r - ro)/we) 



Ai{-C2) 

where Ai is the Airy function and 

Positive constant C2 satisfies 

Ai'(~C2) = 0, 



(27) 

(28) 
(29) 



5.4 Variations in the disc outer radius 

We analyse the eccentricity growth rate due to the 3:1 reso- 
nance as a function of the disc outer radius, or equivalently 
the clearance between the resonance and disc outer edge 
Arc ~ — T-rcs, with all other parameters obtained from the 
standard model. For larger discs, other resonances can play a 
role, such as the 2:1 at r = 0.61 or Arc = 0.15, but we ignore 
such efifects. The solid curve in Fig. [7] plots the eccentricity 
growth rate for the superhump mode as a function of Arc. 
If the disc resides sufficiently far inside the resonant radius, 
that is for Arc < —Wics ~ —0.07, then there is little overlap 
between the disc and the resonant region (i.e., where s(r) is 
substantial; see equation (|19p ). The growth rates under such 
conditions are quite small. The growth rate rises rapidly for 
the disc edge close to, but inside the resonant radius (the 
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Figure 5. Eccentricity distributions for the fundamental free ec- 
centric mode for systems with standard parameters, but different 
values of the dimensionless disc thickness H/r at the disc outer 
edge. The solid curves are determined by numerical solutions to 
equations ((TJ and (I17I I , and the dashed curves are the analytic so- 
lutions for a cool disc given by equation I I27I I. The lowest to high- 
est dashed and solid curve pairs have values of H/r = 0.005, 0.01, 
and 0.02, respectively. 



Arc 

Figure 7. Eccentricity growth rate due to the 3:1 resonance 
in units of Oj, plotted against Arc the difference between the 
disc outer radius and resonance radius. In the standard model 
Arc — 0.034. The solid line is based on smooth continuation of the 
superhump mode. The dashed line is the growth rate for a higher 
order mode. The dotted curve is the 3:1 resonance instability 
factor s(rres -I- Arc) of equation |(TJ in units of Qy^ . Notice that for 
Arc ^ 0.12 this alternative mode dominates, since it has a higher 
growth rate than the superhump mode. 




0.01 - 

0.00 ^ ' 

0.000 0.005 0.010 0.015 0.020 

H/r 

Figure 6. Eccentricity precession rate in units of Q^, plotted 
against H/r, the disc thickness-to-radius ratio at the disc outer 
edge, for the standard model. The dashed line is given by the 
analytic formula, equation I I30I I for a cool disc, the solid line is 
determined by numerical solutions to equation JTJ, and the dot- 
ted line is the precession rate of a free particle at the disc edge, 
rog(ro). 

peak of the dotted line), that is for —0.05 < Arc < 0. Con- 
sequently, the growth rate can be quite sensitive to the disc 
edge location. 

We have idealised the gas distribution to be a smooth 
function of position that is sharply truncated at the disc 
edge. In a more realistic description, the disc density dis- 
tribution tapers near the disc edge. The degree of overlap 
between the disc density distribution can be delicate in the 
case of superhump binaries. The reason is that the disc is 
truncated by tidal forces at a radius that is close to the 3:1 
resonance. Furthermore, the extent of overlap with the reso- 
nant region will likely depend on the disc turbulent viscosity. 
We would qualitatively expect in the superhump binary case 
that the overlap, and consequently the eccentricity growth 
rate, would increase as the disc viscosity increases, since the 



disc tendency to spread outward is increased. Just that type 
of behavior was found in the simulations by Kley el al (2008). 

For the standard model. Arc — 0.034 and the growth 
rate is nearly optimal. For Arc 0.05 the growrth rate drops 
because the peak amplitude of the mode occurs near the disc 
edge and is therefore further away from the resonant region 
(as seen in equation (|19p ). This can be seen in comparing for 
example the solid lines in Fig.|3]and Fig. |S] (that corresponds 
to Arc — 0.14). The former has a stronger overlap with the 
resonance than the latter. So the growth rate in the standard 
disc case (Arc — 0.034) is higher. 

For Arc 0.12, the superhump mode is no longer the 
fastest growing. Such a large disc is unlikely to occur in typ- 
ical superhump binaries, although it could occur in more 
extreme mass ratio systems. Under such a situation, the 
growth rates plotted by the dashed line in Fig. [7] domi- 
nate. The dashed line corresponds to a higher order eccentric 
mode. This mode (actually e^(r), see equation p9|l ) has a 
better overlap with s(r), as seen in Fig.[8l and consequently 
grows faster than the superhump mode. In the absence of a 
companion, this higher order mode would decay faster in the 
presence of viscosity than the superhump mode would, due 
to its smaller scale spatial variations. Other modes could be 
expected to dominate under different conditions. 



5.5 Variations in disc sound speed for 

0.01 <h<0.1 

We consider variations from the standard model in the disc 
gas sound speed or equivalently h, leaving all other param- 
eters fixed. Fig. |9] shows that the eccentricity growth rate 
decreases with the dimensionless disc sound speed, as mea- 
sured by dimensionless disc thickness h. Such behavior was 
found in simulations by Kley et al (2008). As seen in Fig. 1101 
this result can be understood in terms of the effects of the 
sound speed on mode confinement. As the disc sound speed 
increases, the mode spreads further until for h ~ 0.05 the 
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Figure 8. Eccentricity distributions for a system with standard 
parameters, but having disc outer radius To = 0.6 or Arc — 0.14. 
The solid line plots the superhump mode, while the dashed line 
plots a higher order mode, corresponding to the rates plotted in 
solid and dashed curves respectively in Fig. [7] The dotted curve is 
the 3:1 resonance instability factor s(r) of equation Q in units of 
Qb- In this case, the higher order mode overlaps better with the 
instability factor s{r) than the superhump mode and consequently 
has a faster growth rate, as seen in Fig. [T] 



Figure 10. Eccentricity distributions for cases with standard pa- 
rameters, but having different values of dimensionless disc sound 
speed h, the disc thickness-to-radius ratio at the 3:1 resonance. 
The values of h from the lowest to highest solid curves are 0.01, 
0.02, 0.05, and 0.1, respectively. The dotted curves are the 3:1 
resonance instability factors s(r) of equation ((1} in units of Q^, 
for the same set of h values, where the peak values decrease with 
increasing h. 



0.14 




0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 
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Figure 9. Eccentricity growth rate in units of fit, plotted against 
h, the disc thickness-to-radius ratio at the 3:1 resonance, for the 
standard model. The growth rate is higher at lower h because 
the eccentric mode becomes more confined near the resonance, as 
seen in Fig llOl 



eccentricity is fairly uniform. The spreading is not simply 
a consequence of the broadening of the instability factor 
s(r) with h (see equation HI])), since it does not spread as 
rapidly as the mode does with increasing h. Notice that for 
large h, the growth rate in Fig. [9] approaches a constant 
value. This behavior can be understood by the fact that in 
warmer discs the disc eccentricity becomes uniformly dis- 
tributed in radius over the entire disc. We can then crudely 
take We ~ Tres and ercs ~ fimax in equation (I19p and obtain 
a growth rate of 0.02r2b that is in rough agreement with the 
figure for h ~ 0.1. 



5.6 Variations in disc sound speed for h < 0.01 

We analyse the behavior of the eccentricity growth rates and 
eccentricity distributions for cooler discs than were consid- 
ered in Section FS.SI Fig. Illl shows that the growth rate rises 
abruptly with decreasing h, for h < 0.007. Fig. [12] shows 
that the nature of the eccentricity distributions change for 
sufficiently cool discs, since they detach from the disc outer 
edge. In the limit of a very cool disc, the resonance width 
Wrcs ~ h'^^^rrcs becomes smaller than Ar^, the distance from 
the resonance to the disc edge. The distributions are then 
centered on the resonance. This description is valid provided 
that the resonance lies inside the disc, rrcs < ro- The rapid 
growth rates for cool discs are a consequence of the narrow- 
ing of the width of the eccentricity distribution (see equation 

m)- 

If rrcs > ro, the growth rate will drop to zero with 
decreasing sound speed, since there will be little overlap be- 
tween the disc and the resonance. We assume in the remain- 
der of this subsection that the resonance does lie within the 
disc, as follows from the parameters of the standard model. 

We develop an analytic model in the very cool disc limit. 
In that limit, we have that \dE/dr\ » \E\/r and conse- 
quently equation U]) reduces in lowest order to 

in{rros)'^ + 3{x)E{x)^iLoE{x), (32) 

where x = ^j2/'y{r — r^cs ) /H{r res ) • The boundary conditions 
for E{x) are 

£(-oo) = _B'(oo) = (33) 
and the normalization condition is chosen to be 

£■(0) = 1. (34) 
For a cool disc such that H <C w^cs, we approximate 
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the gaussian representation of s{x) given by equation ((9)1 as 



0.8 r 



X 



s{x) 2 

\/27rit;rc 

The solution is given by 



1 ( H 



2\w 



and 



where 



E{r) — exp 



*X 



(35) 



(36) 



2'KWr, 



1 - 



(l + ») 



w. 

WrcB 



— I (37) 



We = 2' 



.5/8^1/4^1/6 



7 



f2 lt)r( 



1/4 



(38) 



The eccentricity e(r-) = \E{r)\ is of gaussian form that 



is centered on the resonance and iAj{r) = (r 



07- 



The real radial velocity determined by equation ^ is 
u'{r, (p) = -rOexp (-(r - T-rcs)^/We) sin ((r - r^^s)'^ /w1 - 4>). 
For r < Trcs, the phase term (r — rrcsY /w1 ~ (j) \a constant 
for increasing r and decreasing 0. The waveform is then 
trailing for r < r^cs and leading for r > Vres- The width uie 
involves a competition between the resonant growth factor 
X and the gas sound speed, through its dependence on H 
and uircs- For Wiaa given by equation (|15p . it follows that 
We/Wrcs goes to zero in the limit that the gas sound speed 
(or H) goes to zero. That is, as the gas sound speed goes to 
zero, eccentricity distribution becomes concentrated into a 
region that is much smaller than the width of the resonance 
uires- From equation (|37|) . we see that eccentricity growth 
rate, —Im{ui), approaches s(rres) for a cool enough disc that 
We <SC Wres. That is, the eccentricity growth rate approaches 
the growth rate of the instability at the resonance, as is 
consistent with equation (|19p . 

Curiously, the detailed properties of the eccentricity dis- 
tribution depend on the form taken for s{x). If instead of a 
gaussian, we adopt a 'top-hat' form, then the results are 
somewhat different. Consider s{x) defined by 



X 



for \/'y/2 \ x\H — \r — rros] < Wrcs and zero otherwise, 
solution to equation (|32|) is then given by 



E{r) 



7r(r 



2Wr, 



for \r — Ticsl < Wias and zero otherwise, and 



(39) 



The 



(40) 



(41) 



Unlike the gaussian case, the width of the eccentricity distri- 
bution does not get smaller than the width of the resonance 
in the very cool disc limit. In this case we have w{r) = 0, 
unlike the gaussian case. While the gaussian case involved 
wrapped modes that were leading (trailing) for r > rrcs 
(r < Tros), the top-hat case contains equal amounts of lead- 
ing and trailing waves at each radius. The growth rate is 
—Im{uj) = x/(2wros) = s(rrcs). So again the eccentricity 
growth rate approaches the growth rate of the instability at 
the resonance. 
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Figure 11. Eccentricity growth rate in units of fit, plotted 
against h, the disc thickness-to-radius ratio at the 3:1 resonance, 
for the standard model. This plot extends Fig. [9] to colder discs. 
For h < 0.007, the growth rate increases more rapidly with de- 
creasing h, as the eccentric mode detaches from the disc outer 
edge and becomes centered on the resonance. 
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Figure 12. Eccentricity distributions, normalised by e(rrcs), for 
cases with standard parameters, but having different values of di- 
mensionless disc sound speed h, the disc thickness-to-radius ratio 
at the 3:1 resonance. The values of h from the lowest to highest 
solid curves at r = 0.5 are 0.005, 0.007, and 0.01, respectively. 
The dotted curves are the 3:1 resonance instability factors s{r) 
of equation JTJ in units of for the same set of h values, where 
the peak values decrease with increasing h. 



5.7 Variations of the inner boundary 

We consider here the influence of both the location of the 
inner boundary and the form of the inner boundary con- 
dition on the eccentricity distribution and growth rate. For 
the cases described above with the standard model or colder 
discs, the inner boundary location typically has little infiu- 
ence. In such cases, the eccentricity distribution smoothly 
approaches zero with decreasing r far from the inner bound- 
ary. But for warmer discs, as may occur in the protostellar 
context, the eccentricity distribution is nonzero close to the 
inner boundary, as seen in Fig. 1101 We consider then the 
case of a standard model disc having h — 0.1. The solid line 
plotted in Fig. [TJ] shows that the growth rate, even in this 
case, is insensitive to the inner boundary radius, provided 
that it is reasonably small, ri < 0.1. 
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Figure 13. Eccentricity growth rate in units of Qj, plotted 
against disc inner radius r\ for the standard model with h = 0.1. 
The solid (dashed) line is for the E = { d^E = 0) inner bound- 
ary condition. The growth rates converge as the disc inner radius 
goes to zero. 



In the protostellar disc context, the disc may not extend 
to the central object, since there may be a hole near the disc 
center. It is then not obvious that the hard wall boundary 
condition (jlip should apply. We consider then the effects of 
applying the inner boundary condition drE = 0, the same 
boundary condition we applied at the outer boundary. To 
explore these effects, we consider an expansion of equation 
ll]) for small r with the standard model. We take r to be di- 
mensionless, normalised by the binary separation. We apply 
equation (|17|) in the small r limit and obtain 



r^E"(r) + 1.5rE'{r) - 0.9E{r) = 0. 



(42) 



The solution is written as 

E{r) = ci{r'" +C2r~''), (43) 

where si ~ 0.731 and S2 ~ 1.23. Constants ci and C2 are 
determined by the inner boundary condition and the re- 
quirement that the solution join onto the solution outside 
of this region of small r. The hard wall (E — 0) boundary 
condition requires C2 = — rf^^"^. In the limit that n goes to 
zero, we have that E oc r"^ . 

The inner boundary condition drE = requires that 
C2 — sir"^^"^ /s2- We expect ci to be of order unity, since 
\E\ ~ 1 where r of order unity. We then have that 

E{rO = C3ri'\ (44) 

for ri small, where cg is a constant of order unity. Therefore, 
we have that the solution with the drE = inner bound- 
ary condition approaches the solution with the hard wall 
E = inner boundary condition in the limit that ri goes to 
zero. That is, both cases have the same E{ri) = value in 
this limit. Consequently, the eccentricity distributions and 
growth rates are the same for both boundary conditions in 
the limit that the inner radius goes to zero. Fig. [13] shows 
that the growth rates are nearly the same for the two bound- 
ary conditions with small ri and that they match as ri goes 
zero. Fig. [14] verifies that equation (I44p holds well for nu- 
merical solutions of eccentric modes when ri is small. 

5.8 Variations of the resonance width 

For the standard model, the resonance width given by equa- 
tion (|15|) evaluates to w^cs — 0.034. Fig. [15] shows how the 
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Figure 14. Eccentricity at the inner disc edge at a function of 
the inner disc radius r; for the standard model with h = O.l. The 
solid line is the result of numerical solutions for equations |[TJ and 
JTtJ with the condition drE = applied at both the inner and 
outer boundaries. The dashed line is the prediction for ri <C 1 
based on equation l|44| l with C3 = 3.82. 
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Figure 15. Eccentricity growth rate in units of fl^, plotted 
against resonance width w with all other parameters given by the 
standard model. The growth rate declines somewhat for larger 
widths, as the resonance instability distribution s(r) spreads be- 
yond the disc outer edge. 



eccentricity growth rate varies with w;res, where all other pa- 
rameters are fixed at their standard values. For small Wics, 
the growth rate approaches a constant value as the resonance 
instability function s(r) approaches the form of a Dirac delta 
function, as was used by GO06. The properties of the eccen- 
tricity distribution are insensitive to the resonance width. 
Recall that for the standard model, the resonance is located 
at distance Arc — 0.034 ~ uircs from the disc outer edge. The 
portion of the resonance distribution s(r) that lies outside 
the disc does not contribute to eccentricity growth. Conse- 
quently, the growth rate declines for w^cs ^ Arc. 



6 SUMMARY AND DISCUSSION 

We have investigated the evolution of discs subject to the 
eccentric instability caused by the 3:1 resonance. The 3:1 
resonance is believed to be responsible for the eccentricity 
of discs in superhump binaries and the excitation mechanism 
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is understood in terms of a mode-coupling model (L91). The 
eccentricity growth rate is affected by the disc eccentricity 
distribution, principally by the radial width of the eccentric 
region and by the eccentricity at the resonance (see equation 
([19])). We explored the disc eccentricity distribution and its 
effect on the growth rate, as a function of system parameters 
using a linear equation developed by GO06, based on a two- 
dimensional disc model. The effects of turbulent viscosity on 
the eccentricity were ignored. 

For typical superhump binary system parameters, the 
disc settles into a mode, the superhump mode, in which the 
eccentricity is concentrated in the outer parts of the disc 
(see Figs. [T]and[3]). The mode results from the excitation 
of fundamental free eccentric mode. For the standard model 
defined in Section [S] the eccentricity growth rate is robust, 
since the mode overlaps well with the resonance and is con- 
centrated near it. If the disc is somewhat smaller or larger 
than for the standard model, then the growth rate of the 
superhump mode is lower (Fig. [7|. Similarly, if tidal forces 
lower the disc density in the region near the resonance, then 
we expect the growth rate to be reduced. In addition, the 
disc tapering needed to avoid Rayleigh instability at the 
outer edge could cause some density reduction in the reso- 
nant region and lower the growth rate. 

The width of the eccentricity distribution is controlled 
by a competition between pressure forces that act to spread 
the mode and the eccentricity growth and differential pre- 
cession that act to confine it (Figs. [3l and llOp . The growth 
rate decreases with the increasing disc sound speed or disc 
thickness (see Fig. [9]) and increases with binary mass ratio 
somewhat faster than quadratically (Fig. Other eccen- 
tric modes can be excited under different conditions, such 
as larger or colder discs (Figs. [8land ll2|) . 

Once the superhump feature is observed in binary sys- 
tems, it is likely that the exponential growth phase, which 
occurs for small eccentricity, has ceased. In that case, the 
disc eccentricity may then be best described by the funda- 
mental free eccentric mode, the free mode excited by the 
superhump instability. Equation (|30p describes the free pre- 
cession rate of a cool disc, subject to the approximations of 
this model. 

Discs in young binary star systems of extreme mass ra- 
tio can similarly undergo eccentric instability. The eccentric- 
ity growth rates in units of the binary frequency are likely 
lower than in the superhump binary case, since protostellar 
discs are considerably warmer as measured by the dime- 
nionless thickness h (see Fig. [9]). Star- planet systems could 
similarly undergo such instability, although it likely involves 
the participation of other resonances that lie closer to the 
planet (Kley & Dirksen 2005, D'Angelo et al 2006). In that 
case also, we expect the eccentricity of the inner disc to be 
fairly broadly distributed because the mass ratio is much 
smaller and the disc thickness is larger than in the super- 
hump binary case. 

The effects of the disc vertical thickness, turbulent vis- 
cosity, and nonlinearity, omitted in this study, should be 
explored in the future (e.g., Ogilvie 2001). 
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